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A dependence of the neutrino masses on the dark energy scalar field could provide a solution to 
the why now problem of dark energy. The dynamics of the resulting cosmological model, growing 
neutrino quintessence, include an attractive force between neutrinos substantially stronger than 
gravity. We present a comprehensive approach towards an understanding of the full cosmological 
evolution including the formation of large-scale neutrino structures. Important effects we account 
for are local variations in the dark energy and the backreaction on the background evolution, as 
well as relativistic neutrino velocities. For this aim, we develop a relativistic A-body treatment of 
the neutrinos combined with an explicit computation of the local quintessence field. At its current 
stage, the simulation method is successful until z ~ 1 and reveals a rich phenomenology. We obtain 
a detailed picture of the formation of large-scale neutrino structures and their influence on the 
evolution of matter, dark energy, and the late-time expansion of the Universe. 



I. INTRODUCTION 

In today's cosmological standard scenario, a cosmolog- 
ical constant A is assumed to explain the observed accel- 
erated expansion of the Universe. Although so far con- 
sistent with all major observational probes [2-[3| > it faces 
two fundamental and unresolved problems. Besides its 
disturbingly tiny value (the cosmological constant prob- 
lem) , it remains miraculous why A has become important 
just recently (the why now or coincidence problem). 

It is thus tempting to think of alternative cosmologi- 
cal scenarios in which these two problems are alleviated. 
In this work, we study growing neutrino quintessence, 
which addresses both of the problems. Growing neu- 
trino quintessence 0, Q relies on two assumptions. First, 
the dark energy component is described by a dynami- 
cal scalar field. Second, the neutrino masses depend on 
this field. Studying the background evolution, it is found 
that the coincidence problem can indeed be solved for a 
neutrino-cosmon coupling somewhat larger than gravity. 

In order to explore the implications of this model, its 
evolution has to be understood also on the perturbation 
level. Linear perturbation theory, however, breaks down 
even at large scales Q • This is due to large overdensities 
in the neutrino fluid becoming non-linear at z < 2 on 
supercluster scales. 

An understanding of the non-linear evolution, al- 
though of utmost importance for confronting the model 
with observational data, is still lacking. In first attempts, 
some aspects of the model have been studied with hydro- 
dynamical and A-body methods [1,11]. Yet, three crucial 
aspects of the model have not been included so far. First, 
as shown by an analytical study of single non-linear neu- 
trino overdensities, local variations of the neutrino mass 
can become very large [Toj . Second, numerical results Q 
show that neutrino particles are likely to reach highly rel- 
ativistic velocities and thereby leave the Newtonian limit. 
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Third, the local mass variations as well as the relativistic 
corrections change the average energy-momentum ten- 
sor of neutrinos; this backreaction effect alters the back- 
ground evolution of the dark energy scalar field. 

It is the purpose of this work to form the basis for 
an adequate A-body based simulation method capa- 
ble of including all major effects of growing neutrino 
quintessence. For the first time, we implement an ex- 
plicit computation of the local dark energy scalar field 
for every time step. This allows us to account for lo- 
cal mass variations and backreaction effects. We present 
a relativistic treatment of the neutrino dynamics as ap- 
propriate for the high velocities of neutrinos occurring 
during the structure formation process. The results of 
linear perturbation theory are used for the initial condi- 
tions. All these aspects are not merely small corrections, 
but decisive for every quantitative prediction. 

The organization of this paper is as follows. We present 
growing neutrino quintessence and derive the equations 
of motion required for our simulation in Sec. [XXJ In 
Sec. IIII1 we discuss strategies how to account for the pe- 
culiar features of growing neutrino quintessence in a nu- 
merical simulation. The numerical results are presented 
in Sees. HVI andlVl At first, we have a look at the forma- 
tion of neutrino structures (Sec. IIV[) . and then explore 
the effects on the evolution of dark energy and matter 
(Sec.|VJ). We summarize in Sec. IVII 



II. GROWING NEUTRINO QUINTESSENCE 

A. Overview 

Growing neutrino quintessence is a possible solution to 
both the cosmological constant and the why now problem 
of dark energy [a, E|. Since it is a quintessence model 
with a dark energy scalar field, the cosmon ip, dark en- 
ergy evolves dynamically. During most of the cosmolog- 
ical evolution, its energy density decays similarly to the 
densities of the other species. The vanishing of the cos- 
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mon potential V(<p) for large values of p> can be rooted 
in the approach to a fixed point with effective dilatation 
symmetry. The scaling solution with dark energy of the 
same order of magnitude as radiation or smaller can then 
explain the tiny overall size of the dark energy density by 
the large age of the Universe. 

In contrast to standard quintessence scenarios, how- 
ever, a coupling between the cosmon <p and the neutrinos 
provides a solution to the coincidence problem as well. 
Due to their small masses, cosmological neutrinos have 
become non-relativistic only recently. This event triggers 
the onset of dark-energy domination in growing neutrino 
quintessence. 

A coupling between the cosmon p and a matter species, 
here assumed to be the neutrinos, is expressed by the 
exchange of energy and momentum. The individual 
energy-momentum tensors do not satisfy a conservation 
equation, only their sum does. Denoting with T a P the 
energy-momentum tensor of the neutrinos and with S a ° 
the energy-momentum tensor of the cosmon field, we 
have 



(1) 



Since no known symmetry requires Q = 0, we generally 
have to expect a non- vanishing coupling. A specific form 
of the coupling proposed by early works [III E3| is 



Q' 



-f3Td a p 



(2) 



with a dimensionless coupling parameter (3 and T 
"!" ,, . using units where the reduced Planck mass Mp = 
(87rG) -1 / 2 is set to unity. Writing T = —pu + Sp^, we see 
that T vanishes as long as the neutrinos are relativistic, 
w v = VvlPv = 1/3- The coupling only becomes effec- 
tive once the neutrinos have become non-relativistic. For 
large values of the coupling (3, the coupling can stop the 
evolution of the cosmon. An almost constant value of the 
cosmon potential leads then to an onset of dark-energy 
domination at recent times, similar to the concordance 
model ACDM. 

On a particle physics level, the coupling, Eq. fl2J), is 
realized as a dependence of the neutrino masses on the 
cosmon field. The coupling parameter /3, in terms of the 
average neutrino mass m„ = m v (p>) ||, now is 
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We further assume that the coupling parameter (3 is con- 
stant, i. e., it does not depend on p. In this case, we have 
the crucial relation 



m v {p) cx e /3(p . 



(4) 



In scenarios with an expansion history similar to ACDM, 
the coupling f} takes large negative values, typically of 
order (3 ~ — 10 2 . The possible couplings to other matter 
species are assumed to be negligible. The self-interaction 



of the cosmon is given by a potential V(ip), for which an 
example is the exponential potential [5], 



V(ip) oc e 



(5) 



where a is a dimensionless model parameter with typical 
values a > 10 to satisfy early dark energy constraints 

003. 

The coupling, strong enough to stop the cosmon evo- 
lution in the background, also crucially modifies the evo- 
lution of perturbations. In fact, if the neutrinos are 
non-relativistic, the cosmon perturbation Sp is approx- 
imately a factor of 2/3 larger than the neutrino-induced 
gravitational potential. The perturbation dp describes 
an attractive force between the neutrinos of order IjFI rj 
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Since realistic scenarios have 



P 3> 1, the extra force leads to a very rapid growth of 
perturbations in the neutrino fluid becoming non-linear 
at z n i ss 1-2, even on large scales Q. 

Analyzing growing neutrino quintessence quantita- 
tively thus requires adequate methods to study the non- 
linear evolution. Perturbations Sp in the quintessence 
field imply local mass variations <5m„ by virtue of Eq. Q • 
These variations can significantly change the averaged 
energy-momentum tensor entering the background equa- 
tions. This backreaction effect of structure formation on 
the background evolution will turn out to be crucial. Fur- 
thermore, the forces during the non-linear evolution are 
strong enough to accelerate the neutrinos to relativistic 
velocities. We thus have to work with a relativistic de- 
scription. For these aims, we shall next derive and collect 
the necessary equations. 



B. Fundamental equations 

In a first step, we will present a fundamental def- 
inition of the model in terms of an action and de- 
rive the basic equations, Eqs. ^ and (0), introduced 
above. In principle, one has to describe all three neutrino 
species together. This can be done in growing neutrino 
quintessence 0). In what follows, however, we adopt the 
simplified working hypothesis of the neutrinos having all 
the same mass. The three flavors of neutrinos enter then 
only in the initial number density of neutrinos. For our 
simulation, we can deal effectively with one species. 

The dynamics of the cosmon field tp and the neutrino 
field ip are described by a usual scalar field Lagrangian C v 
and a Majorana Lagrangian C v respectively, with the pe- 
culiarity that the neutrino mass term m v ipip is assumed 
to depend on tp, 
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C v = --d a p>d a p> - V(<p), 

m v (p))) ip. 



^ 2 



(6) 
(7) 



Here, we use the vierbein formalism to describe the neu- 
trino field in curved spacetime [l5| . The quantities 7" (x) 
are related to the usual Dirac matrices j a (a = 0, 1, 2, 3) 
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by virtue of the vierbein e"(x): 7 Q (x) = -f a e^(x), where 
the vierbein is related to the metric by g a/3 = e"ef^ a6 , 
■q ab = diag(— 1, 1, 1, 1). The action reads 

S = J d 4 x J—g (C + C U + C v ) , (8) 

where Cq includes the remaining cosmological species and 
gravity. A Majorana constraint relates tp to ip. 

We obtain the energy-momentum tensor T a P for the 
neutrino field by the usual definition 

5S V = \J ' d A x^g(Sg^)T^. (9) 

We assume a diagonal metric, which we have, e.g., in 
the Newtonian gauge — otherwise, one employs a gener- 
alization with vierbeins. The energy-momentum tensor 
reads 

T 01 ? = -I^0V a ty + -V (o ^7«V. (10) 

where we have symmetrized in the indices a and /3 (cf., 
e.g., (l5|). The equations of motion following from C v 
are the Dirac equations in curved spacetime, 

7 Q V Q ^ + m^((^)V = 0, 
-VaH a + m u {<p)i> = 0. (11) 

In the absence of the coupling, i. e. in the case m„ = 
const., these equations imply the energy-momentum con- 
servation equation VpT 01 " = for the neutrino field. 

In growing neutrino quintessence, however, the mass 
depends on the spacetime coordinates x via (p(x), and 
the derivative also acts on m„. Inserting the equations 
of motion into VpT a P , we find 

V p T al3 = d a m v (<p) i^V = -d a m v (ip) n v . (12) 

We will give an interpretation of h = —iipip below. Using 
the definition of /3, Eq. ([3]), we may write d a m u (ip) — 
—(3m v (ip)d a ip. In our discussion, an important role is 
played by the trace T of the neutrino energy-momentum 
tensor, 

T = T a a = -m v (ip) h v = — p v + 3p„. (13) 

This eventually shows the form of the energy-momentum 
exchange, cf. Eq. J2]), 

V,gT a/3 = —j3Td a ip. (14) 

A standard computation for the energy-momentum ten- 
sor of the cosmon, 

S af} = d a ipd f) tp + g al3 £ ip , (15) 

leads to 

WpS af3 = +(3Td a ip, (16) 



and one verifies the conservation equation 
V/3 (T a P + S aP ) = for the sum of the neutrino 
and cosmon energy-momentum tensors. 

In the non-relativistic limit, h v introduced above cor- 
responds to the neutrino number density n v . In contrast 
to n u , however, n v does not transform as a scalar. For 
a Lorentz transformation, e.g., n v picks up the volume 
contraction factor I/7. This will become more concrete 
when we represent the neutrino field by effective rela- 
tivistic particles (see Sec. IIII Bj) . 

C. Cosmon dynamics 

The coupling between neutrinos and the cosmon field 
has important impacts on the evolution of both species. 
Variation of the action with respect to (p yields the mod- 
ified Klein-Gordon equation 

V a V a <p-V v (<p) =PT. (17) 

We split into background quantities (spatial averages 
only depending on time) and perturbations (spatially 
varying with vanishing mean), g a p — g a p + 6g a p,(p — 
<p + Sip, and T a @ = T Q/3 + 5T a @ . We assume that the 
metric perturbations Sg a p and the cosmon perturbation 
dip can be treated in linear approximation and that their 
time derivatives are small. 

We choose the conformal Newtonian gauge (see, e.g., 
16]), in which the metric reads 

d.s 2 = a 2 (-(1 + 2*)d?7 2 + (1 - 2$)da: 2 ) , (18) 

with the conformal time 77, the comoving coordinates x, 
the scale-factor a(rj), and the two scalar potentials ^ 
and $. It is now straightforward to evaluate Eq. (jTTJ) . 
Separating into a background part, independent of the 
spatial position, and a linear perturbation part, we obtain 
two evolution equations. 

For the background quantities, we find 

+ 2Ufi + a 2 V v (£) = a 2 /3 [p v - 3p„), (19) 

where primes denote derivatives with respect to r), % = 
a! j a is the conformal Hubble parameter, and we have 
used T = —p v + 3p„. Obviously, the coupling is ineffec- 
tive as long as u>„ = p v /p v = 1/3, i.e. if the neutrinos 
are relativistic. Once the neutrinos turn non-relativistic, 
the right-hand side of Eq. dT9l) stops the further evolu- 
tion of (p. By this mechanism, the model can solve the 
coincidence problem 0, ||| . 

In the perturbations, the equation reads 

&8<p - a 2 V w {Cp)8ip + 2^{(p" + 2H<p') = 

= a 2 f3 5T, (20) 

where the spatial derivatives refer to comoving coordi- 
nates. In the fluid description, ST = —5p v + 35p„. For 
our purpose, however, it is more convenient to calculate 
8T directly from the distribution of particles in our sim- 
ulation. 
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D. Neutrino dynamics 

In this section, we investigate the motion of a neutrino 
particle with a cosmon-depending mass m v {ip). We de- 
scribe a neutrino as a classical, yet relativistic particle 
with world line £ a (r) and four-velocity u a — d^ a /dr, 
where r denotes the particle's proper time, defined via 
dr 2 = — g Q( gd£ a d£^. The energy-momentum tensor of 
this particle is given by 



Comparing this with Eq. (|23|) . we arrive at the equation 
of motion: 



dTm v {tp(£))u a u f} 6 i (x-0, (21) 



where g is the determinant of the metric and 5 A (x) de- 
notes the four-dimensional Dirac delta function. The fac- 
tor l/y/—g ensures the correct normalization of the Dirac 
delta function in curved spacetime by compensating the 
invariant volume form \J—g d 4 x. The one-particle action 
is constructed from the energy-momentum tensor, 



= - / dTm„fo>(0). (22) 



The equations of motion can be obtained by varying S v 
with respect to the particle's path £(t). In the uncoupled 
case, m v = const., this would give the standard geodesic 
equation. The modifications due to the cosmon-neutrino 
coupling are the same that we will find below by using 
energy-momentum conservation. 

In what follows, we will use Eq. (|2 1 [) in the energy- 
momentum conservation equation, Eq. (|14jl. in order to 
derive the equations of motion. Modifications to the stan- 
dard geodesic equation enter in two ways, on the right- 
hand side through the exchange term Q a — —f3Td a tp 
and on the left-hand side through the cosmon-depending 
mass m v (ip). 

Since u a u a = — 1, the right-hand side simply becomes 



du a 
~dr 



T%u"u 



(3d a y + /3u x d x <pu a , (27) 



which describes the deviation from the geodesic motion 
due to the cosmon-neutrino coupling. This equation of 
motion can also be obtained by a conformal transfor- 
mation of the geodesic equation [9|. Let us give a brief 
interpretation of its terms. 

• Tp a .u p u' 7 describes the gravitational effects. In the 
component a = 0, this is the Hubble damping. The 
spatial components give rise to the gravitational 
force. In the Newtonian limit, this would be given 
by V^. In the relativistic case, the situation is 
more complicated involving both potentials, "J and 
$. 

• /3 d a ip is the cosmon-mediated fifth force. In the 
Newtonian limit, this corresponds to an attractive 
force between neutrinos about 2/3 2 times stronger 
than gravity Q. 

• P u*d\<p u a represents a velocity dependent force. 
It can be understood as a consequence of momen- 
tum conservation. Particles are accelerated when 
they move into a direction where they lose mass. 
It also modifies the universal damping term, which 
is no longer proportional to the Hubble parameter 
H, but to (H - /V)- 

It is instructive to consider the case of a static particle, 
u l = 0. Then, the time component of the right-hand side 
of Eq. (|2"7|) vanishes, and the spatial components reduce 
to Pd l p. Thus, a static particle is only affected by the 
spatial variation Vy> and not by ip', and u° = (1 — ^)/a 
does not depend on ip. 



/3Td a ip 



dTm l/ (ip)f3d a ip5 4 (x-£). (23) 



III. STRATEGY AND METHOD 



The left-hand side is 



VpT 01 ? = d p T ali + r^T X13 + T^T"' 3 . (24) 
It is straightforward to show 



dr— KM/)i 4 (!-e). (25) 

OT 



The derivative acting on m v ((p) is dm v {tp)/dT — 
—f3m(p)u x d\ip. Apart from this contribution, i.e. in 
the uncoupled case, Eq. (p4|) reproduces the standard 
geodesic equation. Thus, we find 



v T afS 



&rm v (p) S 4 (x - f) 



dii a 

I + T a pa u p u a - u X d x <p u a 



dr 



(26) 



Although many efforts to describe the structure forma- 
tion process in growin g ne utrino quintessence have been 
made (cf, e.g., 

none of the applied meth- 
ods was capable of drawing a comprehensive picture, 
nor of providing reliable quantitative results once non- 
linearities are strong. 

These challenges motivate the development of a 
new method, specifically designed for growing neutrino 
quintessence. We describe the framework in Sec. IHI Al 
As a fundamental ingredient, we keep track of the per- 
turbed cosmon scalar field allowing us to incorporate lo- 
cal mass differences of the neutrinos. We find a signifi- 
cant neutrino mass suppression once a large fraction of 
the neutrinos is bound in non-linear structures. This in- 
duces a backreaction on the evolution of the cosmological 
background (Sec. IIII Cj) . The initial conditions are drawn 
from the linear results ( Sec. IIII Dj) . We discuss numerical 
issues in Sec. IIII El 
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Framework 



(2) Simulation steps. 



An important goal is to make predictions for the den- 
sity contrasts S v (x) , S m (x) and the peculiar velocity fields 
vf, ec (x), vf£ c (x) of neutrinos and matter respectively. 
These are linked to the gravitational potential and to 
various observables. These fields, however, do not carry 
all the necessary information needed to describe their 
evolution. They merely arise from the first moments of 
the full phase-space distribution functions f u (rj, x l ,Vj), 
fmijIiX 1 ,Vj). In a non-linear evolution, higher moments 
of these distributions become important. It is thus most 
efficient to directly sample the distribution functions by 
a finite number of effective particles, N v and N m respec- 
tively. These particles carry a comoving position x, a 
velocity v = dx/drj, and a rest mass M v (neutrinos), or 
M m (matter). Our method thus bases upon an iV-body 
scheme. 

The motions of these particles depend on the cosmolog- 
ical background, the two gravitational potentials *&(x), 
$(x), and the cosmon ip(x) = <p(r)) + Sp(x). In the New- 
tonian limit and neglecting the local mass variation of 
neutrinos, both gravity and the fifth force can be de- 
scribed as two-body forces. If we want to go beyond 
these approximations, we need to know the explicit val- 
ues of these fields. 

We model the fields VP, $, and Sp as discrete values on a 
three-dimensional grid in a volume V with N c cells. The 
volume is a cube with side length L, every cell has the 
volume (Aa:) 3 = V/N c . We employ periodic boundary 
conditions. For the concrete values chosen, cf. Table |U 



Simulation properties 


Specification 


Box volume V — L 3 


600 3 /i~ 3 Mpc 3 


Number of cells N c 


256 3 


Neutrino particles N v 


2 x 10 7 


Matter particles N m 


2 x 10 7 


Initial redshift 2, 


4 (neutrinos), 49 (matter) 


Final redshift Zf 


1 


Particle properties 


x, v, M„, M m 


Dynamical fields 


*, $, 8tp 


Background quantities 


H, <p, pv,Pv, pm 



TABLE I. The basic parameters and ingredients of the simu- 
lation. 



Let us outline the basic steps of the algorithm we use. 
The details will be given in subsequent sections. 

(1) Initialization of the simulation (cf. Sec. IIII D)l . 

(a) Generate a realization of initial fields <5„, S m , 
t)P ec , «Pf c from the linear spectra. 

(b) Sample the corresponding distribution func- 
tions with N v effective neutrino and N m ef- 
fective matter particles. 



(a) Accelerate and move the effective neutrino 
and matter particles (cf. Sees. IIII Bl and IITTJI) . 

(b) Calculate the potentials ^ , $, and the cosmon 
perturbation Sip (cf. Sec. IIII 51) . Update neu- 
trino masses M v with the new local cosmon 
field if (cf. Sec. lTlTCl) . 

(c) Measure the averages p u , p u . With these 
quantities, evolve the background cosmon (p. 
The Hubble parameter % is evaluated by 
Friedmann's equation 3'H 2 /a 2 = p v + p m + p v 

(cf. Sec. |mc| . 

B. Particles and fields 

As stated above, we assume that the fields 5', $, and 
Sip can be treated linearly and that their time derivatives, 
"J', <f>', and Sip' , are small. This leads to simple algebraic 
equations in Fourier space that allow us to calculate the 
fields from a given distribution of particles. This method 
requires to carry out several Fourier transforms at each 
time step, which can be done efficiently by the use of Fast 
Fourier Transform routines. Let us collect the relevant 
equations. 

The gravitational potential <f>, on subhorizon scales, is 
obtained from the Poisson equation [16j, 



-dp, 



(28) 



where k = \k\ and Sp = Sp u + Sp m + Sp v . 

While we can obtain Sp v from the cosmon field and its 
perturbation, 



Sp v 



ip' Sip 



(29) 



we need to calculate Sp v (and 8p m analogously) from the 
effective particles. 

Since Sp v = — ST°q, we get the contribution of one 
particle at position £ from Eq. (|2T]) . 



1 



1 



/ =jM^S 3 (x 
V9 



0, 



where we have introduced the Lorentz factor, 
_ V^5ooda; _ I 
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dr ^1 - (1 - 2* - 2$)v 2 

and the determinant of the spatial metric, 

g = det{ gij ), v/^ = a 3 (I-3$). 



(30) 



(31) 



(32) 



In the case of matter, we can approximate v 2 <§; 1. In 
contrast, we keep the full relativistic equations for neu- 
trinos since values of v 2 close to one can be reached once 
large neutrino structures form. 
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Finally, the discrete field value Pv{x) at a cell x is 
obtained by summing up the contributions of all the 
particles located inside the cell. Subtracting the mean 
p v yields 5p v {x). The Poisson equation, Eq. ([25]) . can 
now be solved as follows. The sum Sp(x) — 5p v (x) + 
Spm(x) + Sp l p(x) is Fourier transformed and inserted into 
Eq. P5|) . The resulting field $(fc) and its gradient ife<f>(fe) 
are then transformed back to real space, yielding $(x) 
and V$(a;). 

The second gravitational potential ^ equals $ as long 
as anisotropic stress is negligible, i. e. when T l j is small 
for i =/= j. Inspecting the energy-momentum tensor of a 
single particle, cf. Eq. (f2Tj) . the off-diagonal components 
T l j are of order v 2 and thus only important for relativis- 
tic species. In our case, the only contribution comes from 
the neutrinos. Perturbation theory yields the following 
expression for the difference of the two potentials fig : 

k 2 (*-*) = 4 a2 Y.( k -p L -l^)^ 03) 

with the traceless components of the neutrino energy- 
momentum tensor, 

E^T^-^TV (34) 

The field T, l j(x) is obtained from the neutrino particles 
by using the form of the one-particle energy-momentum 
tensor, Eq. (|2T]) . A straightforward calculation yields 

£ % = I^^ 7 M„( l ,„ j -^<s-0. (35) 

The third field, Sip, is calculated from Eq. (Pit))) in 
Fourier space. On the right-hand side, we need 5T = 
ST a a . Again from Eq. ([2T]) . the one-particle contribu- 
tion is 

T=-^^5 3 (x-a (36) 
V9 1 

Given gravitational potentials "J, $, and the cosmon 
perturbation Sip, the acceleration of the effective particles 
is obtained from the corresponding equations of motion. 
In the case of neutrinos, this is Eq. ([27)1 . In the case of 
matter, the right-hand side vanishes, and the equation 
specializes to the usual geodesic equation. 

For matter particles, we make the approximation v 2 <C 
I, and simply find 

^- = -Vy-Mv. (37) 

For the neutrinos, it is adequate to solve the equation 
of motion, Eq. ([27]) . directly for the spatial components 
of the four-velocity u % = 7« l (I — ^)/a. Since this auto- 
matically respects u a u a = —I, the speed of light limit 
is robustly enforced. A differential equation for v itself 
could be unstable in this regard. 



The gradients V^P, V$, 'VSip, occurring in the equa- 
tions of motion are calculated in Fourier space and then 
transformed back to position space like the fields "J, $, 
Sip themselves. This limits the accuracy on scales com- 
parable to the cell size Ax. In particular, forces between 
two particles located in the same cell are completely ne- 
glected. We will investigate the impact of this simplifi- 
cation in Sec. IIII El 

Fourier transformations are also used to estimate the 
power spectrum of perturbation variables, e. g. of the 
neutrino density contrast 5 v (x). The spectrum P v fol- 
lows from the definition, 

(5 v (k)S*(q)) = (2n) 3 P„(k)S 3 (k-q), (38) 

which we discretize on the grid. 



C. Backreaction effects 

The usual way to study cosmological dynamics is to 
perform a complete split between the background and 
the perturbation evolution. First, averaged energy- 
momentum tensors and the unperturbed Friedmann met- 
ric are used to find the evolution of the background quan- 
tities. Second, the evolution of perturbations on the pre- 
calculated background is investigated. 

This procedure neglects the modifications to the back- 
ground evolution due to the presence of perturbations, 
the so-called backreaction. In the standard ACDM 
case, this is a reasonable ap pro ximation since gravita- 
tional backreaction is small 18:]. In growing neutrino 
quintessence, in stark contrast, the backreaction plays a 
decisive role. For coupled quintessence models with cou- 
plings of super-gravitational strength, a significant im- 
pact of the non-linear structure formation on the back- 
ground evolution is expected [H, [1(3, [H, EH ■ 

We now illustrate the importance of backreaction in 
our specific model and explain how we account for it in 
the simulation. The background dynamics of growing 
neutrino quintessence differ from ACDM because they 
include the joint evolution of the cosmon tp and of the 
neutrinos described by their energy-momentum tensor 
T°p. The equation describing this evolution is the mod- 
ified Klein-Gordon equation, Eq. (TT5]). Its right-hand 
side is given by the quantity T = T a a = —p v + 3p„. 

In the standard procedure, neglecting backreaction, 
one would estimate T from the evolution of the averaged 
energy-momentum tensor. Once the neutrinos have be- 
come non-relativistic, the trace is simply given by the 
average energy density, T = —p v , whose evolution fol- 
lows, cf. 0, 

p v + 3Hp v = -fippy. (39) 

This equation tells us that T will depend on the evolution 
of the averaged cosmon tp but does not take into account 
local variations Sip. 
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These local perturbations, however, are very important 
when we regard the evolution of neutrino structures. A 
large part of the neutrinos is concentrated in non-linear 
structures where Sip takes negative values. This leads to 
a systematic suppression of the masses M v oc exp(— /3Stp). 
We shall see that Eq. (|39[) thus leads to a wrong estimate 
of T. 

The correct method to estimate T taking backreac- 
tion into account is to use the exact expression, which 
we obtain from the one-particle contributions, Eq. (|36p . 
Performing the spatial averaging, we get as a sum over 
the particles p, 

J v d 3 x^gT -1 ^ jljyy 

This expression underlines the influence of the pertur- 
bations. The quantity T is suppressed by two essential 
effects. First, we clearly see the dependence on the parti- 
cles' actual masses determined by the local cosmon field. 
Second, once the fifth force has accelerated the particles 
to relativistic velocities, the Lorentz factors 7 > 1 lead 
to an additional suppression. 

One could ask how large the inconsistency becomes if 
one neglected the backreaction effects. For this purpose, 
we run our simulation without modifying the background 
evolution due to the influence of perturbations. Techni- 
cally, we use Eq. (j3"9")l for the background evolution. In 
order to quantify the inconsistency, we measure the av- 
erage, Eq. (|40[) . in every time step. The comparison be- 
tween the two estimates of the same quantity —T = —T a a 
are shown in Fig. [TJ In addition, we show the actual av- 
erage p v of the energy density. 
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-T, actual average 
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FIG. 1. Inconsistency between the assumed background evo- 
lution (green dashed) and the actual average taking into ac- 
count the locally varying cosmon field (blue solid). For illus- 
tration, the figure includes the actual average p„ of the energy 
density (gray dot-dashed). 

For early times, as long as perturbations are small, the 
two estimates agree well. But once non-linear structures 
form, the suppression of —T becomes obvious. This error 



carries over to the evolution of ip by Eq. (fT9)) and thereby 
to the whole subsequent evolution of the background. 

It is thus important to incorporate the backreaction 
effects and to evolve the background and the perturba- 
tions in parallel. In every time step, we thus use the 
actual average, Eq. (|40| , for evolving the cosmon field (p 
bv Eq. (|T9"]). 

We further see the importance of relativistic correc- 
tions. Without these, we would have — T = p u . Figure Q] 
shows that this approximation becomes invalid at late 
times. In order to determine the expansion rate H, we 
measure the actual average p v in our simulation box at 
every time step, instead of using background equations. 
The rate T-L is then obtained from the Friedmann equa- 
tion m 2 = Y, t m 2 - 



D. Initial conditions 

Our simulations use the results from linear perturba- 
tion theory to draw initial conditions. The choice of an 
initial redshift z^ is motivated by two restrictions. First, 
as long as neutrinos are highly relativistic (for z > 5), 
one cannot neglect the higher moments in their phase- 
space distribution function. Second, linear perturbation 
theory becomes invalid at z « 2. We choose Zi — 4 where 
the equation of state w v has fallen to about 10 -2 while 
the neutrino perturbations are still small. 

Although Zi = 4 is a good choice for the neutrinos, 
matter perturbations have to be treated non-linearly 
much earlier. We thus start the TV-body treatment of 
matter at z — 49. During the stage from z — 49 to 
z = 4, matter evolves according to gravity on a growing 
neutrino quintessence background. 

Typically, ./V-body codes for CDM simulations use a 
displacement field to find growing mode initial conditions 
for the particles in the Zel'dovich approximation (see, 
e.g., [lo, HU). In the more complicated case of grow- 
ing neutrino quintessence, however, we prefer to use the 
results of the linear calculation to draw the initial condi- 
tions directly. 



1. Initial random fields 

Before we distribute neutrino particles, we have to 
draw random perturbation fields, namely, the den- 
sity contrast 5 v (r)i,x) and the peculiar velocity field 
v^ ec (r]i,x). We focus on scalar perturbations, i.e., the 
peculiar velocity v? cc is related to a scalar velocity per- 
turbation v s via vf, cc (r],x) = Vw s (?7, x). 

Let us consider generic scalar perturbations 6A(rj,k), 
SB(r), k) in Fourier space. They are solutions of the per- 
turbation evolution equations as given in Q. These are 
linear and only depend on the absolute value k = |fc|. 
Hence, they allow for basis solutions Ak(i]), Bk(v)- We 
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consider the adiabatic mode and write 
(5A( V ,k)\_ (h) fA k ( V ) 



(41) 



The coefficient a(k) characterizes the concrete realiza- 
tion, whereas the basis solutions Ak{r]), Bk(rj) are merely 
attributes of the differential equations. The two-point 
statistics is then 

(8A( V , k)SB*( V , q)) = A k {r,)BM Hk)a*{q)) . (42) 

Assuming that, at an early cosmological time ?7 pr i m , the 
perturbations were described by the primordial spectrum 

-Pprim(^), it follows 



(a(k)a*(q)) = (27r) 3 P prim (fc) <5 3 (fc - q) 



(43) 



We assume Gaussian primordial perturbations with a 
Harrison-Zel'dovich spectrum 



2tt 2 

-Pprim(fc) = ^3 A s 



^pivot 



(44) 



where A s denotes the scalar amplitude, n s the spectral 
index, and fc p i V ot the pivot scale. 

A random realization of the fields SA(r]i,x) and 
8B{rii,x) used as initial fields in the simulation is now 
obtained in three steps. First, a(k) is drawn as a Gaus- 
sian field with the statistics given in Eq. (|4"3")l . Second, the 
linear code [7( is used to obtain the basis solution Ak(r]i) 
and Bk{i]i). Since the linear code works in synchronous 
gauge, a gauge transformation to the Newtonian gauge 
is necessary (cf., e.g., [Hj]). Finally, Eq. (|4T|) determines 
the perturbation fields 5A{rji, k) and 5B(r)i, k) in Fourier 
space, which can be transformed to real space. 

The procedure above also applies to matter starting 
at z = 49. A consistent realization of the two species 
requires to use the same coefficients a(k) for both. 



2. Discrete fields and particle distribution 

In a numerical implementation, real and Fourier space 
are discretized. We will now explain how to draw initial 
conditions on a discrete grid {ki}. Therefore, we first 
have to derive a discretized version of Eq. (|43[) . In par- 
ticular, we work with a discrete Fourier transform (DFT) 
rather than the continuous version. 

First of all, on the right-hand side, the Dirac delta 
function is approximated by the Kronecker delta, 

(2 7 r) 3 P prlm (fc i ) 5 3 {k l - kj) « (2 7 r) 3 P prim (fc l ) 



= P 



prim 



(Afc) 3 
(fc,)L 3 ^, (45) 



where we have used Ak — 2tt/L. On the left-hand side of 
Eq. (|4"3"1) . we replace the quantities a(ki) by the discrete 
Fourier quantities on related to real space by a DFT. 



The relation between a(ki) and on becomes obvious 
when we discretize the Fourier integral. We write fcj = 
*2, ^3) Afc and Xj — (ji, j'3) Ace with Ax = LJn 
and the number n of cells per dimension, i.e. n 3 = N c . 
It follows 



a(ki) — J d 3 ict(x)e" 



TrX>(*i)< 



,-2-n-i {iij-L+i-2j-2+izh)/ n 



N 
V 



(46) 



where we have used the definition of the three- 
dimensional DFT. 

We arrive at the discretized version of Eq. (|4"3")l . 



N 2 



(47) 



expressing the variance of the Gaussian random number 
&i . Only approximately half of the N c numbers a, are in- 
dependent. We have to respect the hermiticity condition 
ensuring that the transformed field a(x) takes on only 
real values. In this way, we obtain our initial random 
fields on the grid. 

Once the fields 8 v (x), S m (x), vf, ec (x), and v^ c (x) are 
calculated, we have to represent them by a distribution 
of iV-body particles. We freely choose a total number 
of A^-body neutrino and matter particles, N u and N m , 
respectively. 

For the neutrinos, we have to account for thermal ve- 
locities v th , which, at z% — 4, are small but not negligi- 
ble. They are described by a Fermi-Dirac distribution 
/ (u th ), cf., e.g., The total particle velocity is then 
v = v th + v pec . We approximate the phase-space distri- 
bution function of neutrinos by 

= ^ fo(v - vl cc (x)) (1 + *„(*)). (48) 

In this way, the averaged particle velocities at each cell 
reproduce the peculiar velocity field. 

The distribution function implies the number density 



n u {x) 



d\f 1/ (x i ,v j ) = ^(l + 5 1/ (x)). (49) 



At each grid point x, we distribute the rounded number 
\n{x) (Air) 3 J of particles and correct the error statisti- 
cally by the addition of further particles. To reduce shot 
noise, particles are put at random positions within the 
pixel volume. The velocity v of a particle in a pixel x 
is drawn from the distribution fo(v — i;P ec (£c)). In order 
to enforce the correct local average, we draw particles 
pairwise with opposite thermal velocities. 

For matter, the procedure is analogous with the only 
difference that thermal velocities are negligible: 



fm (^ j Vj ) 



V 



5 3 (v 



,,pec 



(x)) (1 + 6 m (x)). (50) 
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We have checked this procedure by comparing the 
spectra estimated from the resulting particle distribution 
with those obtained by the linear code. 



E. Numerical issues 

The strategies presented in the previous sections take 
care, in principle, of all relevant effects in growing neu- 
trino quintessence. As an important caveat, the method 
encounters numerical instabilities once very concentrated 
structures have formed. We will discuss this problem in 
the following section. Further, we will address the influ- 
ence of the limited resolution on our results. 



1. Instability for z < 1 



In Sec. IIII B[ we collected several equations that al- 
low us to estimate the potentials \f r , $, and Sp> from the 
distribution of particles. These calculations, however, al- 
ready require knowledge of the potentials, e. g. via the 
metric coefficients, cf. Eq. (132]) . We cannot resolve these 
mutual dependencies exactly. 

As long as the time steps are chosen small enough, 
it is a good approximation to use, in these cases, the 
potentials of the previous step. This is the procedure we 
generally apply. 

This can fail, however, if small errors accumulate in 
such a way that the evolution becomes significantly inac- 
curate or even unphysical. In fact, once the formation of 
neutrino structures is advanced, this problem occurs in 
the estimation of Sp. This limits our capability to simu- 
late growing neutrino quintessence for late cosmological 
times, z < 1. 

In more detail, Eq. (|20[) for Sip contains a source term 
ST = ST a a , which, in turn, is sensitive to the local neu- 
trino mass variations oc exp(— (3Sp>). Formally, Eq. (|20[) 
can be written in the abstract form 



AS<p(x) = f(5p{x);x), 



(51) 



with a highly non-linear function /. 

Attempts to discretize the Laplacian and to apply 
Newton's method to directly solve this non-linear equa- 
tion are numerically intractable due to the large number 
of cells. 

Another approach is to consider the iteration 



AS<p {n+1) (x) = f(Sp {n) (x)-x). 



(52) 



If it converges to a fixed point Sp* , this is clearly a solu- 
tion of Eq. (f5"Tj) . As a starting point S<p(°\ one may use 
the result of the previous time step. 

The convergence of the sequence Sp^ n \ however, is not 
guaranteed. We apply this iterative scheme until z = 1, 
since, for later times, the sequence starts to diverge. 

In order to get a rough idea how the behavior of the 
sequence depends on the physical conditions, we consider 



a very simplified configuration. We assume the idealized 
situation in which all the neutrinos in the volume V are 
concentrated in a dense structure of comoving size / small 
enough to be approximated by a Dirac delta shape. Let 
us fix the origin of the coordinates x at the position of the 
structure. The structure is assigned a total mass AfW cx 
exp(— f3Sp^ (0)) in the n-th iteration. For convenience, 
we take care only of the leading effects, i. e., we drop the 
metric perturbations and neglect V vv in Eq. (l20l) . 
The resulting expression for T reads 



a A 7eff 



(53) 



where 7 c ft arises from the Lorentz factors of the neutrinos 
inside the structure. The iteration prescription, Eq. ([52]) , 
in Fourier space yields, for k ^ 0, 



c (n+i) 1 a 2 



(54) 



Transformed back to real space and evaluated at the 
structure position x = 0, 



2na 1 7eff ' 



(55) 



where we have performed a cut-off of the Fourier in- 
tegral at the scale of the structure size fc max — n/l. 
This value determines the next mass estimate 

M (n+1) 

exp(-^> +1) (0))- 

If we had started with the correct solution Sp>*, the 
mass M* would remain fixed in the iteration. Let us 
assume we had started with a small error, 
M* (1 + £(")). Then we find, at linear order in e^ n \ the 
size of the error in the next step, 



T (n+1) _ 



f3 2 M* 



r (n) 



2wa l-fes 



(56) 



The iteration can only converge if the errors decrease, 
i. e., if the absolute value of the prefactor is smaller than 
unity. 

Of course, the above discussion is only illustrative. 
Nevertheless, it clarifies that the strong coupling /3 and 
large concentrations M*/l of neutrinos in structures are 
working against the iteration. 



2. Resolution 

The equations of motion for the effective neutrino par- 
ticles, Eq. (|27[) . are not reproduced by mere two-body 
forces. Instead, they require the knowledge of the fields 
^(x), $(a;), and Sp(x) on a grid. This differs signifi- 
cantly from pure CDM simulations, which obtain a high 
resolution by using two- body interactions for small-scale 
forces as in GADGET-2 [22j |. Since we do not employ 
an adaptive mesh, our resolution is limited by the con- 
stant size of a grid cell Ax. While this clearly affects the 
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accuracy of our method at small scales (comparable to 
Ax), our approach still seems appropriate for studying 
the model on large scales at which neutrino structures 
first form. 

As a check, we have run our simulation with a lower 
resolution, N c — 128 3 , corresponding to twice a bigger 
cell size Aa; 4/i _1 Mpc. Apart from N c , all other pa- 
rameters as well as the initial random realizations have 
been chosen identical. In Fig. [21 we show the dimension- 
less neutrino spectrum, 

A 2 ,(fc) = ^P,(fc), (57) 

at redshift z = 1 obtained from the two simulation runs. 
As expected, the reduction of N c leads to a loss of power 



cosmological neutrinos are not directly observable, the 
indirect influence of neutrino structures on cosmological 
observables is essential and will be the topic of Sec. [V] 

The model parameters used in our simulation are listed 
in Table [TXJ These are exemplary parameters resulting in 

Model Cosmology 

/3 = -52 Ho = 70 km/s Mpc -1 

a = 10 A s = 2.3 x 10~ 9 

ml = 2.3 eV fc pi vot = 0.05 Mpc" 1 

Q ( l = 0.15 n s = 0.96 
Q.% = 0.60 

TABLE II. Parameters for growing neutrino quintessence and 
primordial perturbations. 
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FIG. 2. The neutrino spectrum A„ at z — 1 for two different 
resolutions. The vertical dashed line marks the scale 1/Ax m 
0.25 /i/Mpc roughly corresponding to the cell size of the low- 
resolution grid. 

on small scales. On the other hand, we find satisfy- 
ing agreement on large scales k < 1/Ax « 0.25 /i/Mpc. 
These results indicate that we can obtain a robust pic- 
ture of large-scale neutrino clustering from our simula- 
tions. When studying the properties of individual neu- 
trino structures in Sec lIVB] the resolution of small scales 
will be more important. We will then again use the low- 
resolution simulation with N c = 128 3 to quantify the 
effect. 



IV. NEUTRINO STRUCTURES 

Using the strategies presented in Sec. IIII1 we are able 
to study growing neutrino quintessence at the non-linear 
level. We first investigate qualitatively how the strong 
fifth force leads to the formation of large neutrino struc- 
tures, Sec. IIV Al Within these structures, the local mass 
variation becomes an important effect as we will see in 
Sec. IIVBI The effects of neutrinos accelerated to rela- 
tivistic velocities are presented in Sec. IIV CI Since the 



a rather large present-day neutrino mass. Very different 
parameter choices are possible, including a time-varying 
coupling constant j3 = /3(<p). The values m°, 0°, and 
would be reached if no backreaction effects were taken 
into account. They do not characterize the state at z — 
of the full cosmological evolution including backreaction. 

A. Growth of structure 

We run the simulation with the cosmological parame- 
ters listed in Table Hi] and the simulation parameters of 
Table HI Recalling that the cosmon-mediated fifth force 
felt by the neutrinos is a factor of about 2/3 2 ps 5 x 10 3 
stronger than gravity, it is expected that large non-linear 
neutrino structures will form. The numerical results we 
get in our simulation agree with this expectation. We 
can follow the precise time evolution of the structure for- 
mation process. We give an overview of the evolution 
from a = 0.25 to a — 0.5 in Fig. G3 The images show 
snapshots of the number density field n v (x) of neutrinos 
in the simulation box (periodic boundary conditions) at 
intervals Aa = 0.05. We will now discuss the main stages 
of neutrino structure formation. 

Until a « 0.3, no large non-linear structures have 
formed. The perturbations can still be described lin- 
early, and our results reproduce the linear calculation 
0- We give an impression of these perturbations and 
their growth by showing a spherical, two-dimensional 
section of the three-dimensional field. Although mainly 
linear, the growth already is very fast. From a — 0.25 
toa = 0.30, the perturbations have grown by a factor of 
about 3 to 4. The particle velocities, however, are still 
far from relativistic. The Newtonian limit would be ap- 
plicable, and we will also see that backreaction effects are 
small. 

From a ~ 0.3 to 0.4, large-scale non-linear inhomo- 
geneities are forming. Most cells in the simulation box 
are already empty of effective neutrino particles, which 
increasingly concentrate in large and dense filament-like 
regions. During this process, the velocities of effective 
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a = 0.25 





FIG. 3. The evolution of non-linear neutrino structures in the simulation box, L = 600 /i _1 Mpc. The lower four figures show 
regions in the simulation box where the number density of neutrinos is a factor of > 5 higher than in the background. We see 
the process of continuing concentration. In the beginning, the structures were still linear. We show, in the upper figures, a 
two-dimensional section of the number density (the color range — blue to red — goes from to 5 times the background value). 



neutrino particles increase and often reach relativistic 
values. The Newtonian limit breaks down, and the fully 
relativistic equation of motion becomes essential. The 
large density inhomogeneities lead to non-negligible lo- 



cal variations Sip of the cosmon and in turn to important 
local mass variations. These modified masses together 
with the Lorentz factors 7 > 1 cause backreaction effects 
according to Sec. IIII CI In fact, we shall find in Sec. IV Al 



12 



that a w 0.4, or z ~ 1.5 is a characteristic time where 
the evolution of the cosmological background has started 
to feel the impact of large inhomogeneities. 

After a w 0.4, approximately spherical neutrino struc- 
tures form (as already predicted in [I3|)j mainly at the 
intersections of the large filament -like non-linearities. In 
the course of the evolution, these structures become even 
increasingly concentrated and spherical. The spherical 
shapes can be regarded as a hint for virialization. In the 
simulation cells, the number densities of effective parti- 
cles exceed the average value by factors up to 10 5 . The 
structures are, however, very different from matter struc- 
tures. They produce local cosmon inhomogeneities Sip 
resulting in locally varying neutrino masses. In the struc- 
tures' centers, the neutrino mass is heavily suppressed, 
cf. Sec. HVBl 

In Sec. IIII E| we explained that our strategy to calcu- 
late the cosmon perturbations Sip becomes problematic 
once the concentration of neutrino structures becomes 
very large. For this reason, we cannot follow the sub- 
sequent cosmological evolution for a > 0.5, i.e. z < 1. 
We expect the overall picture to remain stable, i. e. a 
collection of neutrino structures with very large central 
overdensities. Due to the attractive fifth force, some of 
these structures are expected to merge if they are close 
enough to each other. 



B. Individual structures 

We have found several distinct neutrino structures in 
our simulation box at z = 1. In the following, we have 
a closer look at one of the most pronounced structures. 
In Fig. 21 we show a two-dimensional slice through our 
simulation volume located at the center of the structure. 
Since the structure is almost spherically symmetric, we 
can characterize its shape by a radial profile. We show 
the neutrino number density and the local neutrino mass 
as a function of the physical radial distance r pn from the 
center of the structure. A concentrated core is clearly 
visible with a neutrino number density contrast of 2 x 10 D 
in the central cell. 

A very interesting phenomenon of growing neutrino 
quintessence are the local neutrino mass variations. Since 
the cosmon-mediated fifth forth leads to a substantial 
neutrino clustering in regions of negative cosmon pertur- 
bations Sip, the neutrino mass is typically strongly sup- 
pressed in overdense regions. This suppression is clearly 
visible inside our exemplary structure, cf. Fig. |4] Wc 
find that the neutrino mass at the innermost core of 
the structure is approximately one order of magnitude 
smaller than outside the structure. This result stresses 
the relevance of the local cosmon perturbations. It re- 
mains an open question how the structures evolve for 
z < 1. The idea that the neutrino mass profile inside the 
structures may become constant and independent from 
the evolution of the cosmological background has been 
raised in (Toj . 
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FIG. 4. The uppermost figure shows a slice through the num- 
ber density field of neutrinos (in multiples of the average n„ ) . 
The two lower figures show the number density profile n v (r) 
and the mass profile m v {r) of the structure. The dashed lines 
show the results of a simulation with a lower resolution. The 
shaded regions indicate the size of a grid cell for each resolu- 
tion. 



Our simulation method does not resolve the dynamics 
below the size of a cell, which is around 2 /i _1 Mpc in co- 
moving units. Hence, at small distances from the center, 
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the quantitative results are affected by rather large er- 
rors. In order to estimate the uncertainties, we show the 
corresponding results obtained from a simulation with a 
twice as large cell size (dashed lines). The number den- 
sity profile shows a similar shape but is significantly less 
concentrated. Still, the mass suppression reaches similar 
values. 

We have seen that the neutrino structures are still be- 
coming increasingly concentrated in the final stage of our 
evolution (cf. Fig. [3]) . This is most visible for the largest 
structures and less significant for smaller ones. As an ex- 
ample, we show how the number density profile of one of 
the smaller, but still very pronounced structures evolves 
from a = 0.45 to a — 0.5 in Fig. [5j We have verified that 
the structure is not undergoing merging processes in the 
considered time range. Apart from a moderate transfer 



average over the contributions of single particles p to the 
density p v = -T° , Eq. flU, yields 
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FIG. 5. The figure shows the number density profile n v of an 
isolated neutrino structure at a = 0.50 (blue solid), a — 0.47 
(green dashed), and a = 0.45 (gray dot-dashed) as a function 
of the physical distance from its center. 

of neutrinos from the outer regions of the structure to 
its inner core, no significant changes in the profile are 
visible. 



C. Relativistic effects 

The motion of the effective neutrino particles in our 
simulation respects the relativistic laws of motion, cf. 
Sec. Ill D[ which remain valid for velocities approaching 
the speed of light. This is essential since the neutrinos feel 
a strong acceleration due to the fifth force. We observe 
in our simulations that the Newtonian laws of motion 
assuming small velocities become inadequate as soon as 
non-linear clustering begins, z ss 2 to 1.5, cf. Sec. IIV Al 

The importance of the relativistic treatment becomes 
manifest in the evolution of the average equation of state 
Wv = Pu/Pu- With the help of the methods explained 
in Sec. MI B[ we can calculate the average pressure p v 
and energy density p v of the neutrino fluid from the dis- 
tribution of particles in our simulation. Performing the 



Pv 



J v d 3 xVgT°o _ 1 
Jv 



d^7r=^E^(^))- (58) 



The average pressure p v then follows from the calculation 
of f = f a a by Eq. (@D1), p v = (p v + f )/3. 

The evolution of the equation of state w v is shown in 
Fig. |B1 The figure shows a steep increase of w v , grow- 
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FIG. 6. The increase of the equation of state io„ due to rel- 
ativistic velocities in the structure formation process. The 
dashed horizontal line marks the limit w = 1/3 for highly 
relativistic particles. 

ing several orders of magnitude between z = 2.5 and 
z = 1, exceeding the value of w v » 0.1. The velocities 
of many neutrino particles in the box are close to the 
speed of light at this stage of the evolution. The fig- 
ure also shows the evolution of w v obtained by solving 
the homogeneous background equations, i. e. neglecting 
backreaction effects. If neutrinos were not accelerated by 
the cosmon-mediated fifth force, their equation of state 
would continue to decrease. The pronounced oscillatory 
features are attributed to the neutrino mass variation; 
the cosmon field oscillates around the minimum of its 
effective potential (cf. @). 

We observe a striking discrepancy between the back- 
ground calculation and the evolution of the actual av- 
erage. The consequences for the evolution of the back- 
ground cosmon field (p will be discussed in Sec. IV A[ 

The presence of relativistic neutrinos is also interesting 
with regard to gravity. Our tools allow us to calculate 
the anisotropic stress induced by neutrinos, which is the 
source for the difference between the two gravitational 
potentials "J/ and Sec. IIIIB1 We find that the effect 
is most pronounced in the vicinity of neutrino structures 
and becomes negligible on very large scales. In Fig. [71 
we show the field $ - $ at z = 1 on the same two- 
dimensional slice already used in Fig. [4] In the regions of 
neutrino overdensities, the field $ — \& shows anisotropic 
patterns with amplitudes of the order 10 -5 . This is, com- 
pared to $ and $ themselves, a 1% to 10% effect. 
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FIG. 7. The difference $ — 9 of the two gravitational poten- 
tials, scaled by a factor of 10 , at z = 1. 



V. IMPACT ON DARK ENERGY AND 
MATTER 

Once the dynamics of growing neutrino quintessence 
can be described reliably, a confrontation with observa- 
tional data is the next step. At the current stage, our 
simulation method allows us to follow the cosmological 
evolution only until z ft; 1, and we have not explored the 
parameter space of growing neutrino quintessence. It is 
thus not yet possible to provide constraints on the model 
parameters. 

Nonetheless, it is insightful to discuss some remarkable 
effects linked to observations. This regards the expansion 
history of the Universe (depending on the perturbation 
evolution due to strong backreaction) , Sec. IV A| and the 
evolution of matter perturbations, Sec. IV Bl 



A. Backreaction and quintessence 



own right. In growing neutrino quintessence, the evolu- 
tion of the dark energy is intimately connected, and hence 
sensitive, to the evolution of the neutrinos, cf. Sec. [TTJ In 
particular, the cosmological event of the neutrinos be- 
coming non-relativistic serves as a trigger stopping the 
further evolution of the cosmon and leading to an epoch 
of accelerated expansion. We have seen, however, that 
the fifth force accelerates neutrinos again to relativistic 
velocities in the course of the non-linear evolution. This 
reduces the neutrinos' capability of stopping the cosmon 
evolution. The onset of accelerated expansion is thus ex- 
pected to shift to later times. 

More precisely, the source term in the background evo- 
lution of the cosmon in Eq. (flUf is proportional to the 
trace T. This trace was shown to be a sum over particle 
contributions cx m^/j, Eq. (j4*U|). So, both the effect of 
relativistic particles, 7 > 1, and the mass suppression go 
into the same direction. 

For a quantitative investigation, we calculate the de- 
celeration parameter q = -a" a/a' +1. For the expan- 
sion of the Universe to accelerate, q has to take neg- 
ative values. We show the evolution of q measured in 
our simulation, i.e., including backreaction according to 
Eq. (|4"U|). compared with the expectation based on the 
a-priori averaged background equations, cf. Eq. p9[) . in 
Fig. [5J We observe that the deceleration parameter, al- 
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We argued in Sec. lIII Cl that two essential effects in the 
structure formation process lead to an important backre- 
action of the perturbations on the background evolution. 
These are the systematic suppression of neutrino masses 
within structures, m v (ip(x)) < m v ((p), and the relativis- 
tic Lorentz factors, 7 > 1. Numerically, we have seen 
that these effects are indeed strong. The neutrino mass 
inside structures is suppressed by an order of magnitude, 
Fig. [31 and the neutrino equation of state w v grows con- 
siderably at late times, Fig. [5] Consequently, the true 
averaged energy-momentum tensor T a p of the neutrino 
fluid significantly differs from the idealized calculation 
based on a homogeneous fluid of non-relativistic neutri- 
nos. 

The strong backreaction is not only interesting in its 



FIG. 8. Evolution of the deceleration parameter q in our 
simulation — including the non-linear backreaction effects — 
(blue solid), compared to the result obtained by the back- 
ground equations (green dashed). 

ready very close to zero without backreaction, is indeed 
still at a relatively large value q ~ 0.2 at z = 1. Although 
our simulation ends at z = 1, we expect the phase of 
accelerated expansion to start later due to backreaction 
effects. This is because we can expect that the two main 
contributions, the mass suppression and the relativistic 
Lorentz factors, remain important in the subsequent evo- 
lution. 

By means of Friedmann's equations, we can express 
the deceleration parameter q as a sum over contributions 
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of the different species i, 

q ~ m 2 



(59) 



The cosmon contribution to the deceleration thus is 



■(p v + 3p v ) oc ip' 2 - a 2 V(ip), 



(60) 



which is a comparison of the kinetic energy and the po- 
tential energy (the potential is normalized to V(0) — 
0.9 x 10~ 120 Mf, with the reduced Planck mass M P ). We 
can use V(ip) = V((p) in linear approximation in Sp. The 
difference between the actual average V(ip) and V((p) is 
way below the percent level. 

Due to Eq. ([60|). we can understand the evolution of the 
deceleration parameter, Fig. [FJ by examining the back- 
ground evolution of the cosmon. We show the equation 
of state w v as well as the average values tp and (p 1 in 
Fig. [9] 

In the lowermost plot, we see that once backreaction 
effects become important (at redshifts around z w 1.5), 
the time derivative (p 1 takes significantly larger values 
than it would be expected neglecting backreaction. This 
shows that, in fact, the neutrino fluid is less effective in 
stopping the cosmon evolution. This is also reflected in 
the plot of tp, which continues to grow although, neglect- 
ing backreaction, a mere oscillation around a very slowly 
increasing value has set in at z < 1.5. As a consequence, 
the equation of state w v is, due to the backreaction, fur- 
ther away from the cosmological constant value w = — 1. 

We conclude that the study of non-linear structure 
formation in growing neutrino quintessence turns out to 
be crucial for the understanding of its correct background 
evolution as well. A realistic model with the correct frac- 
tion of dark energy today might need a readjustment of 
the model parameters. 



B. Matter perturbations 



2 



-0.5 



-0.6 



10 



-10 



with backreaction 
no backreaction 




2 2.5 

redshift z 





with backreaction 


/ 1 \ 
/ ' \ 


no backreaction 


Jf- 1 \ 
1 \ 
\ 1 \ 

•' \ 

- 1 / \ 
\ ' \ 
1 / \ 
\ / \ 









1.5 



2.5 

redshift . 



3.5 



The dynamics of matter perturbations in our simula- 
tion are only described by the gravitational force oc 
and the Hubble damping. The effective matter parti- 
cles are accelerated according to the Newtonian law of 
motion, Eq. (|37[) . and we do not differentiate between 
baryons and dark matter. Matter is not coupled to the 
neutrinos or to the cosmon except via gravity. 

The gravitational effects due to the neutrinos and the 
cosmon include their gravitational potential and their im- 
pact on the expansion history. In particular, the large 
neutrino structures investigated in Sec. IIVI induce gra- 
dients in the gravitational potential leading to an ad- 
ditional acceleration of matter particles. We thus ex- 
pect that the amplitudes of matter perturbations are en- 
hanced compared to the ACDM case [1, [23[ • 



FIG. 9. Evolution of the equation of state of quintessence, 
Wip, the background cosmon field <p, and its time derivative 
(p 1 (including backreaction effects: blue solid; compared to the 
background equations: green dashed). 

In order to study this effect, we follow the evolution of 
the dimensionless matter spectrum, 

The concrete values of A m (k) of course depend on the 
model parameters we have chosen for growing neutrino 
quintessence, Table [TT1 Hence, in order to isolate the im- 
pact of non-linear neutrino structures on the growth of 
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matter perturbations, we normalize A m (fc) by the case 
where matter only feels its own gravitational potential. 
This means, we have started another simulation run fol- 
lowing the evolution of only matter particles in a universe 
with the expansion history of unperturbed growing neu- 
trino quintessence. Since the expansion history is, for 
the chosen set of model parameters, similar to ACDM, 
we somewhat imprecisely label the corresponding ampli- 
tudes by A^ CDM (/e). 

The quotients A m (fc)/A^ CDM (fc) at different redshifts 
are shown in Fig. [TUJ 



1.14 




0.01 0.1 1 
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^Era 2 - (63) 

i 

We show the evolution of Ui for two large scales in 
Fig.DU 
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FIG. 11. The enhancement of matter bulk velocities 
j^^acdm in subboxes of vo i um e / 3 for I = 37.5 /i _1 Mpc 
(blue solid) and I = 75fa~ 1 Mpc (green dashed). 



FIG. 10. The evolution of the relative matter spectrum 

A / a ACDM 



At early times, z > 2, the distribution of the neutrinos 
is still very smooth and matter fluctuations grow only due 
to their own gravitational potential. At later times, we 
observe an enhancement of matter fluctuations at large 
scales, where the neutrino fluid is clustering, cf. Fig. [3] At 
z = 1, the amplification exceeds 10% on the largest scales 
of our simulation, whereas it is still below the percent 
level at small scales. 

By means of the continuity equation, the change S' m 
of the matter density perturbation is linked to the veloc- 
ity field v^ c . If the density perturbations show a steep 
increase, as we observe it on large scales, large peculiar 
velocities must be present. 

We can measure the matter bulk flow in a subvol- 
ume V of our simulation defined as the average peculiar 
velocity, 



f v d 3 xy/gv% 
Jv d 3 x y/g 



(62) 



Next, we introduce the variance of the bulk flow Uf on 
different comoving scales I as a measure of the expected 
bulk velocity in volumes of size I 3 . Technically, we divide 
the simulation box into n subvolumes V% of equal size I 3 = 
L 3 /n. We then obtain the average peculiar velocity t>£^ 
within the boxes i and compute the variance according 



Indeed, the enhancement of bulk flows is much more 
pronounced than the increase of density perturbations. 
This feature of growing neutrino quintessence has been 
studied in [H, [23J]. The increase of bulk flows reaches 
factors of about 1.5 to 2 at z = 1, and it is growing. 

The large bulk flows are the consequence of large-scale 
gradients in the gravitational potential Whereas mat- 
ter perturbations are still linear on scales k < 0.1 h/Mpc, 
the neutrino density field is highly inhomogeneous on 
these and larger scales. It comes as no surprise that the 
non-linear neutrino structures dominate the large-scale 
gravitational potential. In order to quantify this, we es- 
timate the power spectrum P^{z, k) of the gravitational 
potential. The dimensionless spectrum 



* 2 (z,fc)^ ^P*(z,fc) 



(64) 



is a measure for the fluctuation variance of the cosmolog- 
ical gravitational potential on spatial volumes rs tt 3 /k 3 . 
We compare the large-scale evolution of *&(z, k) (includ- 
ing neutrino structures) with 1 i /ACI3M (z,k) (neglecting 
neutrino structure formation) in Fig. 1121 Since $ = \& on 
large scales to a good approximation, cf. Sec. IIV C[ the 
corresponding plot for <fr(z, k) is visually identical. 

Once the non-linear effects of neutrino structure for- 
mation become important at z < 2, the large-scale gravi- 
tational potential grows drastically compared to ACDM. 
Since the matter perturbations have only grown mod- 
erately under the influence of neutrino structures until 
z = 1, see Fig. [TUl this effect is due to the neutrino- 
induced gravitational potential. Consequently, Fig. Q2] 
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FIG. 12. The evolution of the gravitational potential 
*(z,fe)/* ACDM (z,fc) for modes k = 0.05 h/Mpc (blue solid) 
and k = 0.02/i/Mpc (green dashed). 



shows that the large-scale potential of the neutrino struc- 
tures exceeds that of matter by about an order of mag- 
nitude. The large-scale gravitational potential ^(fc), 
Eq. at z = 1 reaches values between 10 -5 and 10 -4 
in our simulation. 

Whereas the large-scale gravitational potentials in 
ACDM are constant during matter domination and very 
slowly decay thereafter, growing neutrino quintessence 
shows a completely different behavior. During the for- 
mation of neutrino structures, the large-scale gravita- 
tional potentials become deeper. This is expected to 
have a drastic impact on the integrated Sachs- Wolfe ef- 
fect, which is sensitive to the time derivative (<F + $)'. 
In this way, the evolution of the gravitational potentials 
might be decisive for scrutinizing the model. A continu- 
ing growth of the gravitational potentials for z < 1 would 
clearly be in conflict with observational results [24| . The 
evolution for z < 1, however, could be very different. 
The virialization of the neutrino structures might stop 
the growth of the neutrino-induced gravitational poten- 
tial. 



VI. SUMMARY 

The cosmological dynamics of growing neutrino 
quintessence show a much higher complexity than the 
standard ACDM cosmology or models of uncoupled dark 
energy. This work has made a crucial step towards a 
consistent simulation of growing neutrino quintessence 
including all its relevant effects. 

The standard methods of linear perturbation theory 
and (Newtonian) iV-body simulations, while very suc- 
cessful in the ACDM case, rely on approximations in 
conflict with the effects of growing neutrino quintessence. 
Linear perturbation theory breaks down on all scales Q , 
and in Newtonian A^-body simulations, the particles' ve- 



locities exceed the speed of light [9j. Furthermore, the 
local variation of the average neutrino mass due to per- 
turbations in the dark energy scalar field could so far only 
be studied in idealized configurations [Toj]. None of the 
applied methods has shed light on the expected backreac- 
tion fl9l ] of the highly non-linear neutrino perturbations 
on the evolution of the cosmological background. 

Wc have thus decided to develop a new method coping 
with these challenges right from the start. This effort 
has led to an A-body based simulation, designed from 
scratch and adjusted to growing neutrino quintessence. 
In particular, we incorporate local variations in the dark 
energy scalar field and respect relativistic dynamics for 
the neutrino component. We run the simulation for an 
exemplary choice of model parameters. Our results show 
that all the effects that had to be neglected in previous 
approaches are indeed important. 

The simulation confirms the formation of large-scale 
neutrino structures (already predicted in [17j), and we 
give a detailed description of their evolution (Sec. lIV Al . 
As a main result, we observe a significant neutrino mass 
suppression of an order of magnitude inside concentrated 
structures (Sec. I IV Bp agreeing with [lfj. This is the 
consequence of significant inhomogeneities in the dark 
energy scalar field, the cosmon (p. Once these become 
too large, our method encounters numerical difficulties 
(Sec. lIIIE|) . which is why we do not follow the cosmolog- 
ical evolution for redshifts z < 1. 

We find that the velocities of neutrinos are accelerated 
to relativistic values during the process of structure for- 
mation. The average neutrino equation of state reaches 
values w v « 0.1 at z « 1 (Sec. IIVC[) . Further, we have 
shown that the relativistic motion induces a difference 
between the two gravitational potentials. In the vicinity 
of large neutrino structures, the difference $ — \& typically 
amounts to ~ 10~ 5 . 

Local mass variations and relativistic velocities of the 
neutrinos lead to strong backreaction effects (Sec. 11110)1 , 
We have demonstrated that the backreaction can no- 
tably modify the late-time expansion of the Universe 
(Sec. IV A[) . It is likely that the onset of the accelerated 
expansion occurs later as compared to a homogeneous 
approximation. 

The evolution of matter perturbations is affected 
by the neutrino-induced gravitational potential, which 
dominates on large scales (Sec. IV B|) . As a consequence, 
large-scale bulk flows of matter are amplified by a factor 
close to two at z = 1. The density field, however, is 
much less affected. On large scales, we have found an 
effect of about 10%. 

The current stage of the simulation method presented 
in this work is already promising, and the efforts to ex- 
tend its range of applicability will continue. When the 
cosmological evolution can be followed until z = for a 
collection of different model parameters, a confrontation 
of growing neutrino quintessence with observational data 
becomes possible. 
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So far, we are only able to identify several effects that 
might be interesting regarding observations. The large- 
scale neutrino-induced gravitational potentials could be 
observable directly via gravitational lensing or indirectly 
via the resulting large-scale bulk flows of matter and the 
enhanced density power spectrum. The time evolution 
of the gravitational potential, on the other hand, leaves 
imprints in the integrated Sachs-Wolfe effect. In con- 
trast to ACDM, we have observed increasing large-scale 
potentials in growing neutrino quintessence during the 
formation of neutrino structures. 

The main motivation to search for alternatives to the 
cosmological constant scenario are the difficulties to 
understand the tiny value of A and the why now problem. 
Any competing model should avoid them in the first 
place. Growing neutrino quintessence offers a mechanism 



to solve the why now problem of dark energy. With a 
rich phenomenology beyond the standard scenario, the 
prospects are promising that the model eventually can 
be put to stringent tests once its non-linear evolution is 
completely understood. 
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